Objective: To forecast a real time series using ETS and ARIMA models.
set.seed(32296622)
myseries <- aus_retail |>
filter(!(`Series ID` %in% c("A3349561R","A3349883F","A3349499L","A3349902A",
"A3349588R","A3349763L","A3349372C","A3349450X",
"A3349679W","A3349378T","A3349767W","A3349451A"))) |>
filter(`Series ID` == sample(`Series ID`,1))
The data set used is a subset of the aus_retail data which contains one numeric variable, i.e. Turnover in $Million AUD. aus_retail is a part of tsibbledata package and can also be loaded with the package fpp3 which contains total 9 packages including tsibbledata.
The data is a time series of class tsibble and the source for the same is Australian Bureau of Statistics, catalogue number 8501.0, table 11.(Robjhyndman, n.d.).
The sample data set contains the Turnover data for the ‘hardware, building and garden supplies retailing’ for the state of ‘South Australia’ for the years ‘1982 to 2018’ (calculated each month from April 1982 to December 2018) and contains 441 observations and 5 columns.
The dataset contain the following variables and table 1.1 displays the variable type as well.
data_dict <- data.frame(
S.No. = c("1", "2", "3", "4", "5"),
Variables = c ("State", "Industry", "Series ID", "Month", "Turnover" ),
DataType = c("Character", "Character", "Character", "Month", "Numeric"))
knitr::kable ( data_dict ,caption = "Data Dictonary") |>
kable_styling(latex_options = c("striped", "hold_position")) |>
kable_paper("hover", full_width = T) |>
scroll_box(width = "100%", height = "300px")
| S.No. | Variables | DataType |
|---|---|---|
| 1 | State | Character |
| 2 | Industry | Character |
| 3 | Series ID | Character |
| 4 | Month | Month |
| 5 | Turnover | Numeric |
myseries
## # A tsibble: 441 x 5 [1M]
## # Key: State, Industry [1]
## State Industry `Series ID` Month Turnover
## <chr> <chr> <chr> <mth> <dbl>
## 1 South Australia Hardware, building and garden … A3349501L 1982 Apr 8.6
## 2 South Australia Hardware, building and garden … A3349501L 1982 May 9.5
## 3 South Australia Hardware, building and garden … A3349501L 1982 Jun 8.4
## 4 South Australia Hardware, building and garden … A3349501L 1982 Jul 10.3
## 5 South Australia Hardware, building and garden … A3349501L 1982 Aug 10.6
## 6 South Australia Hardware, building and garden … A3349501L 1982 Sep 11.5
## 7 South Australia Hardware, building and garden … A3349501L 1982 Oct 10.8
## 8 South Australia Hardware, building and garden … A3349501L 1982 Nov 13.1
## 9 South Australia Hardware, building and garden … A3349501L 1982 Dec 25.4
## 10 South Australia Hardware, building and garden … A3349501L 1983 Jan 9.2
## # ℹ 431 more rows
Table 1.2 displays an yearly average turnover in $Million AUD.
year_average <- myseries |>
mutate(year=year(Month)) |>
index_by(year) |>
summarise(Total_turnover= mean(Turnover))
knitr::kable ( year_average ,caption = "Yearly Turnover Summary") |>
kable_styling(latex_options = c("striped", "hold_position")) |>
kable_paper("hover", full_width = T) |>
scroll_box(width = "100%", height = "300px")
| year | Total_turnover |
|---|---|
| 1982 | 12.02222 |
| 1983 | 11.17500 |
| 1984 | 13.39167 |
| 1985 | 14.46667 |
| 1986 | 12.91667 |
| 1987 | 14.03333 |
| 1988 | 16.21667 |
| 1989 | 20.43333 |
| 1990 | 23.45000 |
| 1991 | 23.94167 |
| 1992 | 28.28333 |
| 1993 | 27.84167 |
| 1994 | 28.45833 |
| 1995 | 27.71667 |
| 1996 | 31.39167 |
| 1997 | 34.41667 |
| 1998 | 39.23333 |
| 1999 | 42.93333 |
| 2000 | 55.35833 |
| 2001 | 70.86667 |
| 2002 | 74.66667 |
| 2003 | 70.90833 |
| 2004 | 69.04167 |
| 2005 | 62.55833 |
| 2006 | 67.86667 |
| 2007 | 75.57500 |
| 2008 | 96.25000 |
| 2009 | 88.49167 |
| 2010 | 85.64167 |
| 2011 | 85.88333 |
| 2012 | 70.88333 |
| 2013 | 66.82500 |
| 2014 | 66.57500 |
| 2015 | 79.94167 |
| 2016 | 82.10833 |
| 2017 | 92.55833 |
| 2018 | 98.35000 |
myseries |>
features(Turnover, list(mean = mean)) |>
arrange(mean)
## # A tibble: 1 × 3
## State Industry mean
## <chr> <chr> <dbl>
## 1 South Australia Hardware, building and garden supplies retailing 51.1
Here the data has been divided into 4 equal section, each containing 25% of the data.
myseries |>
features(Turnover, quantile)
## # A tibble: 1 × 7
## State Industry `0%` `25%` `50%` `75%` `100%`
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 South Australia Hardware, building and garden … 8.4 24.6 53.5 76.7 124.
Figure 1.1 shows that both the ends of this Q-Q plot deviates, and hence, has a fat tail but its center also does not follow a typically straight line. This means that the plot has many extreme values and the values vary throughout as well.
ggplot(myseries, aes(sample = Turnover)) +
geom_qq() +
geom_qq_line(color = "#6C03A0")+
ylab("Sample Turnover Values ($Million AUD)")+
xlab("Theoretical Values")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 1.1: Quantile-Quantile Plot
Figure 1.2 is a histogram and is somewhat similar to Bi-modal data, i.e., it represents 2 peaks. Here, it seems that both peaks have similar densities and the curve is sinusoidal.
ggplot(myseries, aes(x=Turnover)) +
geom_histogram(aes(y=after_stat(density)), fill = "#6EBA6B")+
geom_density(color = "#6C03A0") +
ylab("Density")+
xlab("Turnover ($Million AUD)")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 1.2: Distribution of Turnover
Figure 1.3 shows that:
Trend: There is somewhat an increasing trend visible in the data set.
Seasonality: There is a multiplicative seasonality as there is a change in height of the values with time. There are strong peaks observed at the end of each year.
myseries |>
autoplot(Turnover, color= "#C63F2C")+
xlab("Year/Month")+
ylab("Turnover ($Million AUD)")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 1.3: Time Series
According to figure 1.4,the highest turnovers are observed in December in each year and February has vales towards the lower end. This is possible due to the fact that hardware, building and garden work are outdoor tasks and is preferred to be done during summers and summer rises in December in Australia and starts to end by February.
myseries |>
gg_season(Turnover,
labels = "both")+
ylab("Turnover ($Million AUD)")+
scale_colour_viridis_c()+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 1.4: Seasonal Plot
In figure 1.5 the blue line shows the mean value for each month for all years. This plot confirms our observation from figure 1.4 that the highest turnover was observed in December and lowest in February.
myseries |>
gg_subseries(Turnover, color= "#156A06")+
ylab("Turnover ($Million AUD)")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 1.5: Subseries Plot
Below it is seen that a high value is present for the seasonal features.
myseries |>
features(Turnover, feat_stl)
## # A tibble: 1 × 11
## State Industry trend_strength seasonal_strength_year seasonal_peak_year
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 South Austr… Hardwar… 0.991 0.872 9
## # ℹ 6 more variables: seasonal_trough_year <dbl>, spikiness <dbl>,
## # linearity <dbl>, curvature <dbl>, stl_e_acf1 <dbl>, stl_e_acf10 <dbl>
Here in figure 1.6, the colours indicate different months on the vertical axis. Most of the lag plots have shown positive correlation and hence, confirming strong seasonality in the data.
myseries |>
gg_lag(Turnover, lags=1:24, geom='point', size=0.5)+
facet_wrap(~ .lag, ncol=6)+
ylab("Turnover ($Million AUD)")+
xlab("Lag")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 1.6: Lag plot
Autocorrelation is shown in figure 1.7, i.e. a linear relationship between lagged values of time series. It shows that:
There are peaks observed at every 12th lag, i.e. at the end of each month and hence, there is seasonality.This is also suggested by the fact that an equivalent amount of dip is observed (from the previous lag) at equal intervals.
The *Scalloped” shape observed is also due to seasonality.
This data has a trend as autocorrelations for the small lags are large and positive (as observations close in time are also close in value) and these positive values slowing decreases as lags increase.
myseries |>
ACF(Turnover, lag_max = 48) |>
autoplot()+
ylab("ACF")+
xlab("Lag [1 Month]")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 1.7: ACF
Below is a summary of seven autocorrelation features (in this order): first autocorrelation feature of the original data, sum of squares of first 10 coefficients from original data, first autocorrelation feature from differenced data, sum of squares of first 10 coefficients from differenced data, first autocorrelation feature from twice differenced data, sum of squares of first 10 coefficients from twice differenced data and autocorrelation coefficient at first seasonal lag.
myseries |>
features(Turnover, feat_acf)
## # A tibble: 1 × 9
## State Industry acf1 acf10 diff1_acf1 diff1_acf10 diff2_acf1 diff2_acf10
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 South Aust… Hardwar… 0.957 7.84 -0.143 0.0904 -0.501 0.267
## # ℹ 1 more variable: season_acf1 <dbl>
Transformations and difference is done to make the data stationary. There are 3 steps for this process:
Check the changing variance in the data, if visible, transform the data with an appropriate lambda value. Evaluate if the data is stationary and if no, move to the next step.
Seasonal difference to be carried out to remove seasonality, if the data is seasonal. This step also removes trend sometimes. Evaluate if the data is stationary and if no, move to the next step.
Perform regular difference to remove any trend and anything else leftover.
Figure 2.1 shows a time series by transforming turnover to log(turnover) and it seems that there is still some trend and seasonality in the data, so this is not an appropriate transformation.
myseries |>
autoplot(log(Turnover))+
ylab("Log Turnover")+
xlab("Year Month")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 2.1: Log Transformation
Therefore, performing the Guerrero test ahead to find out lambda value for box-cox transformation.
lambda<- myseries |>
features(Turnover, features = guerrero)
lambda <- pull(lambda)
lambda
## [1] 0.4662629
In figure 2.2 box-cox transformation is performed on turnover with lambda= 0.4662629 and it still seems to have trend and seasonal pattern. The data is not stationary.
myseries |>
autoplot(box_cox(Turnover, lambda))+
ylab("Transformed Turnover (lambda=0.467)")+
xlab("Year Month")+
theme_minimal()+
theme( plot.background = element_rect(fill = "#FFF8F7"),
panel.grid.major.y = element_blank(),
axis.text.x = element_text(color="#993333", angle=45, size=8),
axis.text.y = element_text(colour= "#08007F", size=8))
Figure 2.2: Box-Cox Transformation
Figure 2.3 shows the time plot and ACF together and adds a new plot to the report, i.e. PACF (will be used in section 3).
myseries |>
gg_tsdisplay(box_cox(Turnover, lambda), plot_type = "partial")
Figure 2.3: Timeplot, ACF, PACF with Transform
There are still features of some trend and seasonality visible in the ACF, i.e. data is not stationary and therefore, requires seasonal difference.
Confirming conclusion from figure 2.3 using unit root KPSS test which defines Null Hypothesis Ho= Data is Stationary.
myseries |>
features(box_cox(Turnover, lambda), unitroot_kpss)
## # A tibble: 1 × 4
## State Industry kpss_stat kpss_pvalue
## <chr> <chr> <dbl> <dbl>
## 1 South Australia Hardware, building and garden supplies … 6.91 0.01
Here the p value is 0.01, which actually means that the value is < 0.01 (due to a limit applied) If p < 0.05, we need to difference the data. Hence, Ho is rejected.
Confirming if a seasonal difference is needed.
myseries |>
features(box_cox(Turnover, lambda), unitroot_nsdiffs)
## # A tibble: 1 × 3
## State Industry nsdiffs
## <chr> <chr> <int>
## 1 South Australia Hardware, building and garden supplies retailing 1
It is shown that 1 seasonal difference is needed.
In figure 2.3, ACF shows that, the there is seasonality and is a monthly lag data, and hence, the lag is used as 12 for the first seasonal difference.
myseries |>
gg_tsdisplay(difference(box_cox(Turnover, lambda), lag=12), plot_type = "partial")
<img src=“project_files/figure-html/fig11, fig.cap(”After Transforming and Seasonal Differencing”)-1.png” width=“80%” />
From ?? can be said stationary as it is observed that:
In the ACF, the data falls to 0 quickly after the 11th lag.
The time plot shows unrelated values, i.e. not varying with time.
There seems to be no seasonality or trend left in the data.
This is confirmed by unit root KPSS tests below:
myseries |>
features(difference(box_cox(Turnover, lambda), lag=12), unitroot_kpss)
## # A tibble: 1 × 4
## State Industry kpss_stat kpss_pvalue
## <chr> <chr> <dbl> <dbl>
## 1 South Australia Hardware, building and garden supplies … 0.155 0.1
Here the p value is 0.1 (or > 0.1) and can be said stationary.
Further checking if more differencing is required:
myseries |>
features(difference(box_cox(Turnover, lambda), lag=12), unitroot_ndiffs)
## # A tibble: 1 × 3
## State Industry ndiffs
## <chr> <chr> <int>
## 1 South Australia Hardware, building and garden supplies retailing 0
This shows that no more difference is needed, i.e. no need to perform step 3 of evaluating a regular difference. Hence, the data is now stationary.
ETS
myseries |>
autoplot(box_cox(Turnover, lambda))
view <- myseries
This has trend and seasonality
test <- myseries |>
filter(Month > yearmonth("2009 Dec"))
train <- myseries |>
filter(Month <= yearmonth("2009 Dec"))
Selecting ETS
This is a multiplicative seasonal series. This is also the case for error and it follows a general additive trend. The appropriate models could be (MAM), (MNM) or (MAdM)
fit2 <- train |>
model(MAM= ETS(box_cox(Turnover, lambda)~ error("M")+ trend("A")+ season("M")),
MNM= ETS(box_cox(Turnover, lambda)~ error("M")+ trend("N")+ season("M")),
MAdM= ETS(box_cox(Turnover, lambda)~ error("M")+ trend("Ad")+ season("M")),
ets_auto= ETS(box_cox(Turnover, lambda))
)
report(fit2)
## # A tibble: 4 × 11
## State Industry .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 South Au… Hardwar… MAM 0.00253 -682. 1397. 1399. 1462. 0.191 0.255 0.0370
## 2 South Au… Hardwar… MNM 0.00265 -687. 1405. 1406. 1462. 0.190 0.257 0.0379
## 3 South Au… Hardwar… MAdM 0.00266 -687. 1410. 1412. 1479. 0.189 0.258 0.0376
## 4 South Au… Hardwar… ets_a… 0.167 -662. 1354. 1355. 1411. 0.160 0.227 0.305
glance(fit2)
## # A tibble: 4 × 11
## State Industry .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 South Au… Hardwar… MAM 0.00253 -682. 1397. 1399. 1462. 0.191 0.255 0.0370
## 2 South Au… Hardwar… MNM 0.00265 -687. 1405. 1406. 1462. 0.190 0.257 0.0379
## 3 South Au… Hardwar… MAdM 0.00266 -687. 1410. 1412. 1479. 0.189 0.258 0.0376
## 4 South Au… Hardwar… ets_a… 0.167 -662. 1354. 1355. 1411. 0.160 0.227 0.305
tidy(fit2)
## # A tibble: 65 × 5
## State Industry .model term estimate
## <chr> <chr> <chr> <chr> <dbl>
## 1 South Australia Hardware, building and garden supplies… MAM alpha 0.493
## 2 South Australia Hardware, building and garden supplies… MAM beta 0.000100
## 3 South Australia Hardware, building and garden supplies… MAM gamma 0.329
## 4 South Australia Hardware, building and garden supplies… MAM l[0] 4.12
## 5 South Australia Hardware, building and garden supplies… MAM b[0] 0.0410
## 6 South Australia Hardware, building and garden supplies… MAM s[0] 0.905
## 7 South Australia Hardware, building and garden supplies… MAM s[-1] 0.893
## 8 South Australia Hardware, building and garden supplies… MAM s[-2] 0.969
## 9 South Australia Hardware, building and garden supplies… MAM s[-3] 1.58
## 10 South Australia Hardware, building and garden supplies… MAM s[-4] 1.08
## # ℹ 55 more rows
accuracy(fit2)
## # A tibble: 4 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Sout… Hardwar… MAM Trai… -0.264 3.63 2.48 -0.815 6.28 0.442 0.444 0.160
## 2 Sout… Hardwar… MNM Trai… 0.266 3.59 2.46 0.761 6.25 0.439 0.439 0.160
## 3 Sout… Hardwar… MAdM Trai… 0.218 3.58 2.45 0.521 6.23 0.437 0.439 0.137
## 4 Sout… Hardwar… ets_a… Trai… 0.327 3.23 2.24 0.602 5.90 0.400 0.396 0.0800
fit2 |>
forecast(h=24) |>
autoplot(train, level = NULL)
fit2 |>
forecast(h=24) |>
accuracy(myseries)
## # A tibble: 4 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 MAM Sout… Hardwar… Test -4.54 7.53 6.60 -5.58 7.82 1.18 0.922 0.412
## 2 MAdM Sout… Hardwar… Test 1.33 5.59 4.57 1.15 5.24 0.814 0.684 0.309
## 3 MNM Sout… Hardwar… Test 1.19 5.50 4.50 0.984 5.17 0.801 0.673 0.304
## 4 ets_auto Sout… Hardwar… Test 3.91 6.32 4.84 4.27 5.40 0.862 0.774 0.254
ARIMA
train |>
gg_tsdisplay(difference(box_cox(Turnover, lambda), lag=12), plot_type = "partial")
Finding (p,d,q) and (P,D,Q) d=0, D=1 For q we look at ACF plot and significant lags before 12 is only 1st, can be 2nd too. So q=1 or can be q=2. In PACF, sig lag is on 13th , so p=1 .
Hence (p,d,q)= (1,0,1) (2,0,1) (0,0,1) (2,0,0) (1,0,2) (2,0,2) (0,0,2) (2,0,0)
For seasonal values, we only look at seasonal lags, i.e. multiples of 12. Here we did 1 difference so D=1 looking at ACF, Q=0 as ACF quickly decays to 0. from PACF, only lag at 12 and 24 seem significant so P=2 Hence (P,D,Q)= (2,1,0)
Looking at the time plot, the average is not the center so a constant can be present.
Trying combinations
fit <- train |>
model(
arima101210= ARIMA(box_cox(Turnover, lambda) ~ pdq(1,0,1) +PDQ(2,1,0)),
arima201210= ARIMA(box_cox(Turnover, lambda) ~ pdq(2,0,1) +PDQ(2,1,0)),
arima_auto= ARIMA(box_cox(Turnover, lambda) , trace = TRUE, stepwise = FALSE)
)
## ARIMA(0,0,0)(0,1,0)[12]+c 2383.780219
## ARIMA(1,0,0)(0,1,0)[12]+c 2065.762027
## ARIMA(2,0,0)(0,1,0)[12]+c 2059.718262
## ARIMA(3,0,0)(0,1,0)[12]+c 2061.696775
## ARIMA(4,0,0)(0,1,0)[12]+c 2063.796158
## ARIMA(5,0,0)(0,1,0)[12]+c 2066.853129
## ARIMA(0,0,1)(0,1,0)[12]+c 2213.080292
## ARIMA(1,0,1)(0,1,0)[12]+c 2057.582008
## ARIMA(2,0,1)(0,1,0)[12]+c 2060.456079
## ARIMA(3,0,1)(0,1,0)[12]+c 2063.494362
## ARIMA(4,0,1)(0,1,0)[12]+c 2065.886516
## ARIMA(5,0,1)(0,1,0)[12]+c Inf
## ARIMA(0,0,2)(0,1,0)[12]+c 2140.384920
## ARIMA(1,0,2)(0,1,0)[12]+c 2059.460906
## ARIMA(2,0,2)(0,1,0)[12]+c 2062.645203
## ARIMA(3,0,2)(0,1,0)[12]+c Inf
## ARIMA(4,0,2)(0,1,0)[12]+c 2053.515803
## ARIMA(0,0,3)(0,1,0)[12]+c 2114.322638
## ARIMA(1,0,3)(0,1,0)[12]+c 2061.537588
## ARIMA(2,0,3)(0,1,0)[12]+c Inf
## ARIMA(3,0,3)(0,1,0)[12]+c Inf
## ARIMA(0,0,4)(0,1,0)[12]+c 2104.330096
## ARIMA(1,0,4)(0,1,0)[12]+c 2061.840268
## ARIMA(2,0,4)(0,1,0)[12]+c Inf
## ARIMA(0,0,5)(0,1,0)[12]+c 2066.749516
## ARIMA(1,0,5)(0,1,0)[12]+c 2067.738778
## ARIMA(0,0,0)(1,1,0)[12]+c 2388.912382
## ARIMA(1,0,0)(1,1,0)[12]+c 1988.871122
## ARIMA(2,0,0)(1,1,0)[12]+c 1979.787526
## ARIMA(3,0,0)(1,1,0)[12]+c 1980.346459
## ARIMA(4,0,0)(1,1,0)[12]+c 1980.131926
## ARIMA(5,0,0)(1,1,0)[12]+c 1982.852326
## ARIMA(0,0,1)(1,1,0)[12]+c 2198.536697
## ARIMA(1,0,1)(1,1,0)[12]+c 1977.527624
## ARIMA(2,0,1)(1,1,0)[12]+c 1979.159286
## ARIMA(3,0,1)(1,1,0)[12]+c 1981.058057
## ARIMA(4,0,1)(1,1,0)[12]+c 1982.219936
## ARIMA(0,0,2)(1,1,0)[12]+c 2107.674480
## ARIMA(1,0,2)(1,1,0)[12]+c 1978.521179
## ARIMA(2,0,2)(1,1,0)[12]+c 1980.742735
## ARIMA(3,0,2)(1,1,0)[12]+c Inf
## ARIMA(0,0,3)(1,1,0)[12]+c 2072.693151
## ARIMA(1,0,3)(1,1,0)[12]+c 1980.566226
## ARIMA(2,0,3)(1,1,0)[12]+c 1982.395783
## ARIMA(0,0,4)(1,1,0)[12]+c 2055.766572
## ARIMA(1,0,4)(1,1,0)[12]+c 1980.434290
## ARIMA(0,0,5)(1,1,0)[12]+c 2020.513023
## ARIMA(0,0,0)(2,1,0)[12]+c 2402.184994
## ARIMA(1,0,0)(2,1,0)[12]+c 1980.379405
## ARIMA(2,0,0)(2,1,0)[12]+c 1969.126339
## ARIMA(3,0,0)(2,1,0)[12]+c 1970.287296
## ARIMA(4,0,0)(2,1,0)[12]+c 1971.497779
## ARIMA(0,0,1)(2,1,0)[12]+c 2209.318389
## ARIMA(1,0,1)(2,1,0)[12]+c 1965.551612
## ARIMA(2,0,1)(2,1,0)[12]+c 1968.278252
## ARIMA(3,0,1)(2,1,0)[12]+c 1971.385328
## ARIMA(0,0,2)(2,1,0)[12]+c 2116.577744
## ARIMA(1,0,2)(2,1,0)[12]+c 1967.176326
## ARIMA(2,0,2)(2,1,0)[12]+c Inf
## ARIMA(0,0,3)(2,1,0)[12]+c 2078.503433
## ARIMA(1,0,3)(2,1,0)[12]+c 1969.241921
## ARIMA(0,0,4)(2,1,0)[12]+c 2058.024486
## ARIMA(0,0,0)(0,1,1)[12]+c 2380.639874
## ARIMA(1,0,0)(0,1,1)[12]+c 1942.334373
## ARIMA(2,0,0)(0,1,1)[12]+c 1929.930584
## ARIMA(3,0,0)(0,1,1)[12]+c 1928.128341
## ARIMA(4,0,0)(0,1,1)[12]+c 1927.943266
## ARIMA(5,0,0)(0,1,1)[12]+c 1930.693891
## ARIMA(0,0,1)(0,1,1)[12]+c 2188.074954
## ARIMA(1,0,1)(0,1,1)[12]+c 1924.086536
## ARIMA(2,0,1)(0,1,1)[12]+c 1925.923538
## ARIMA(3,0,1)(0,1,1)[12]+c 1929.181886
## ARIMA(4,0,1)(0,1,1)[12]+c 1929.870287
## ARIMA(0,0,2)(0,1,1)[12]+c 2092.922018
## ARIMA(1,0,2)(0,1,1)[12]+c 1924.641434
## ARIMA(2,0,2)(0,1,1)[12]+c 1927.896355
## ARIMA(3,0,2)(0,1,1)[12]+c 1927.772661
## ARIMA(0,0,3)(0,1,1)[12]+c 2048.064739
## ARIMA(1,0,3)(0,1,1)[12]+c 1926.656888
## ARIMA(2,0,3)(0,1,1)[12]+c 1929.796161
## ARIMA(0,0,4)(0,1,1)[12]+c 2025.849955
## ARIMA(1,0,4)(0,1,1)[12]+c 1927.589656
## ARIMA(0,0,5)(0,1,1)[12]+c 1996.580382
## ARIMA(0,0,0)(1,1,1)[12]+c 2389.348326
## ARIMA(1,0,0)(1,1,1)[12]+c 1946.439456
## ARIMA(2,0,0)(1,1,1)[12]+c 1935.629449
## ARIMA(3,0,0)(1,1,1)[12]+c 1936.081051
## ARIMA(4,0,0)(1,1,1)[12]+c 1934.404582
## ARIMA(0,0,1)(1,1,1)[12]+c 2168.886184
## ARIMA(1,0,1)(1,1,1)[12]+c 1928.914405
## ARIMA(2,0,1)(1,1,1)[12]+c 1935.088220
## ARIMA(3,0,1)(1,1,1)[12]+c 1935.712184
## ARIMA(0,0,2)(1,1,1)[12]+c 2079.958572
## ARIMA(1,0,2)(1,1,1)[12]+c 1930.005122
## ARIMA(2,0,2)(1,1,1)[12]+c 1936.311494
## ARIMA(0,0,3)(1,1,1)[12]+c 2042.263345
## ARIMA(1,0,3)(1,1,1)[12]+c 1931.677145
## ARIMA(0,0,4)(1,1,1)[12]+c 2024.472131
## ARIMA(0,0,0)(2,1,1)[12]+c Inf
## ARIMA(1,0,0)(2,1,1)[12]+c 1943.395224
## ARIMA(2,0,0)(2,1,1)[12]+c 1926.177511
## ARIMA(3,0,0)(2,1,1)[12]+c 1925.136462
## ARIMA(0,0,1)(2,1,1)[12]+c 2179.908932
## ARIMA(1,0,1)(2,1,1)[12]+c 1924.907915
## ARIMA(2,0,1)(2,1,1)[12]+c 1923.472617
## ARIMA(0,0,2)(2,1,1)[12]+c 2093.331345
## ARIMA(1,0,2)(2,1,1)[12]+c 1925.890525
## ARIMA(0,0,3)(2,1,1)[12]+c 2044.325935
## ARIMA(0,0,0)(0,1,2)[12]+c 2368.879469
## ARIMA(1,0,0)(0,1,2)[12]+c 1944.392023
## ARIMA(2,0,0)(0,1,2)[12]+c 1931.983933
## ARIMA(3,0,0)(0,1,2)[12]+c 1930.177234
## ARIMA(4,0,0)(0,1,2)[12]+c 1929.958695
## ARIMA(0,0,1)(0,1,2)[12]+c 2173.021754
## ARIMA(1,0,1)(0,1,2)[12]+c 1926.123889
## ARIMA(2,0,1)(0,1,2)[12]+c 1927.940249
## ARIMA(3,0,1)(0,1,2)[12]+c 1931.217251
## ARIMA(0,0,2)(0,1,2)[12]+c 2084.946350
## ARIMA(1,0,2)(0,1,2)[12]+c 1926.663945
## ARIMA(2,0,2)(0,1,2)[12]+c 1929.975449
## ARIMA(0,0,3)(0,1,2)[12]+c 2041.541538
## ARIMA(1,0,3)(0,1,2)[12]+c 1928.697730
## ARIMA(0,0,4)(0,1,2)[12]+c 2023.172967
## ARIMA(0,0,0)(1,1,2)[12]+c Inf
## ARIMA(1,0,0)(1,1,2)[12]+c 1942.631072
## ARIMA(2,0,0)(1,1,2)[12]+c 1930.332835
## ARIMA(3,0,0)(1,1,2)[12]+c 1930.412742
## ARIMA(0,0,1)(1,1,2)[12]+c Inf
## ARIMA(1,0,1)(1,1,2)[12]+c 1924.467953
## ARIMA(2,0,1)(1,1,2)[12]+c 1929.275549
## ARIMA(0,0,2)(1,1,2)[12]+c 2080.917511
## ARIMA(1,0,2)(1,1,2)[12]+c 1925.565573
## ARIMA(0,0,3)(1,1,2)[12]+c 2040.632120
## ARIMA(0,0,0)(2,1,2)[12]+c Inf
## ARIMA(1,0,0)(2,1,2)[12]+c 1943.946678
## ARIMA(2,0,0)(2,1,2)[12]+c 1928.277684
## ARIMA(0,0,1)(2,1,2)[12]+c Inf
## ARIMA(1,0,1)(2,1,2)[12]+c 1926.402310
## ARIMA(0,0,2)(2,1,2)[12]+c 2084.021278
## ARIMA(0,0,0)(0,1,0)[12] 2445.856124
## ARIMA(1,0,0)(0,1,0)[12] 2068.956068
## ARIMA(2,0,0)(0,1,0)[12] 2061.003309
## ARIMA(3,0,0)(0,1,0)[12] 2062.402596
## ARIMA(4,0,0)(0,1,0)[12] 2064.106671
## ARIMA(5,0,0)(0,1,0)[12] 2067.167357
## ARIMA(0,0,1)(0,1,0)[12] 2253.590114
## ARIMA(1,0,1)(0,1,0)[12] 2058.295754
## ARIMA(2,0,1)(0,1,0)[12] 2060.925053
## ARIMA(3,0,1)(0,1,0)[12] 2063.925288
## ARIMA(4,0,1)(0,1,0)[12] 2066.178110
## ARIMA(5,0,1)(0,1,0)[12] Inf
## ARIMA(0,0,2)(0,1,0)[12] 2167.850728
## ARIMA(1,0,2)(0,1,0)[12] 2059.963226
## ARIMA(2,0,2)(0,1,0)[12] Inf
## ARIMA(3,0,2)(0,1,0)[12] Inf
## ARIMA(4,0,2)(0,1,0)[12] 2054.840529
## ARIMA(0,0,3)(0,1,0)[12] 2134.701085
## ARIMA(1,0,3)(0,1,0)[12] 2062.017429
## ARIMA(2,0,3)(0,1,0)[12] Inf
## ARIMA(3,0,3)(0,1,0)[12] 2051.251800
## ARIMA(0,0,4)(0,1,0)[12] 2120.899119
## ARIMA(1,0,4)(0,1,0)[12] 2062.565653
## ARIMA(2,0,4)(0,1,0)[12] Inf
## ARIMA(0,0,5)(0,1,0)[12] 2075.744698
## ARIMA(1,0,5)(0,1,0)[12] 2060.951421
## ARIMA(0,0,0)(1,1,0)[12] 2453.856057
## ARIMA(1,0,0)(1,1,0)[12] 1993.276870
## ARIMA(2,0,0)(1,1,0)[12] 1981.604540
## ARIMA(3,0,0)(1,1,0)[12] 1981.394696
## ARIMA(4,0,0)(1,1,0)[12] 1980.656113
## ARIMA(5,0,0)(1,1,0)[12] 1983.165778
## ARIMA(0,0,1)(1,1,0)[12] 2260.260611
## ARIMA(1,0,1)(1,1,0)[12] 1978.585515
## ARIMA(2,0,1)(1,1,0)[12] 1979.682551
## ARIMA(3,0,1)(1,1,0)[12] 1981.588279
## ARIMA(4,0,1)(1,1,0)[12] 1982.746398
## ARIMA(0,0,2)(1,1,0)[12] 2158.666976
## ARIMA(1,0,2)(1,1,0)[12] 1979.131187
## ARIMA(2,0,2)(1,1,0)[12] 1981.237863
## ARIMA(3,0,2)(1,1,0)[12] Inf
## ARIMA(0,0,3)(1,1,0)[12] 2114.174306
## ARIMA(1,0,3)(1,1,0)[12] 1981.099805
## ARIMA(2,0,3)(1,1,0)[12] 1982.753323
## ARIMA(0,0,4)(1,1,0)[12] 2089.767674
## ARIMA(1,0,4)(1,1,0)[12] 1981.343797
## ARIMA(0,0,5)(1,1,0)[12] 2049.548873
## ARIMA(0,0,0)(2,1,0)[12] 2454.289229
## ARIMA(1,0,0)(2,1,0)[12] 1984.759783
## ARIMA(2,0,0)(2,1,0)[12] 1970.952047
## ARIMA(3,0,0)(2,1,0)[12] 1971.434295
## ARIMA(4,0,0)(2,1,0)[12] 1972.051039
## ARIMA(0,0,1)(2,1,0)[12] 2260.524224
## ARIMA(1,0,1)(2,1,0)[12] 1966.505474
## ARIMA(2,0,1)(2,1,0)[12] 1968.927371
## ARIMA(3,0,1)(2,1,0)[12] 1972.046684
## ARIMA(0,0,2)(2,1,0)[12] 2165.817813
## ARIMA(1,0,2)(2,1,0)[12] 1967.856077
## ARIMA(2,0,2)(2,1,0)[12] 1971.425185
## ARIMA(0,0,3)(2,1,0)[12] 2123.496834
## ARIMA(1,0,3)(2,1,0)[12] 1969.850724
## ARIMA(0,0,4)(2,1,0)[12] 2097.412076
## ARIMA(0,0,0)(0,1,1)[12] 2445.134786
## ARIMA(1,0,0)(0,1,1)[12] 1952.518763
## ARIMA(2,0,0)(0,1,1)[12] 1935.450460
## ARIMA(3,0,0)(0,1,1)[12] 1931.697734
## ARIMA(4,0,0)(0,1,1)[12] 1930.882179
## ARIMA(5,0,0)(0,1,1)[12] Inf
## ARIMA(0,0,1)(0,1,1)[12] 2253.336877
## ARIMA(1,0,1)(0,1,1)[12] 1927.517053
## ARIMA(2,0,1)(0,1,1)[12] 1928.485450
## ARIMA(3,0,1)(0,1,1)[12] 1931.681633
## ARIMA(4,0,1)(0,1,1)[12] 1932.611365
## ARIMA(0,0,2)(0,1,1)[12] 2157.858148
## ARIMA(1,0,2)(0,1,1)[12] 1927.357305
## ARIMA(2,0,2)(0,1,1)[12] 1930.307521
## ARIMA(3,0,2)(0,1,1)[12] 1930.246600
## ARIMA(0,0,3)(0,1,1)[12] 2113.867591
## ARIMA(1,0,3)(0,1,1)[12] 1929.427834
## ARIMA(2,0,3)(0,1,1)[12] 1932.593064
## ARIMA(0,0,4)(0,1,1)[12] 2087.863580
## ARIMA(1,0,4)(0,1,1)[12] 1930.744523
## ARIMA(0,0,5)(0,1,1)[12] 2045.809209
## ARIMA(0,0,0)(1,1,1)[12] Inf
## ARIMA(1,0,0)(1,1,1)[12] 1955.984521
## ARIMA(2,0,0)(1,1,1)[12] 1939.838234
## ARIMA(3,0,0)(1,1,1)[12] 1939.346458
## ARIMA(4,0,0)(1,1,1)[12] 1936.999786
## ARIMA(0,0,1)(1,1,1)[12] 2258.506511
## ARIMA(1,0,1)(1,1,1)[12] 1932.126805
## ARIMA(2,0,1)(1,1,1)[12] 1936.916024
## ARIMA(3,0,1)(1,1,1)[12] 1937.996621
## ARIMA(0,0,2)(1,1,1)[12] 2159.359073
## ARIMA(1,0,2)(1,1,1)[12] 1932.573012
## ARIMA(2,0,2)(1,1,1)[12] 1938.261547
## ARIMA(0,0,3)(1,1,1)[12] 2116.163758
## ARIMA(1,0,3)(1,1,1)[12] Inf
## ARIMA(0,0,4)(1,1,1)[12] 2091.112530
## ARIMA(0,0,0)(2,1,1)[12] 2455.077637
## ARIMA(1,0,0)(2,1,1)[12] 1952.648104
## ARIMA(2,0,0)(2,1,1)[12] 1932.365742
## ARIMA(3,0,0)(2,1,1)[12] Inf
## ARIMA(0,0,1)(2,1,1)[12] 2262.363875
## ARIMA(1,0,1)(2,1,1)[12] Inf
## ARIMA(2,0,1)(2,1,1)[12] Inf
## ARIMA(0,0,2)(2,1,1)[12] 2167.581875
## ARIMA(1,0,2)(2,1,1)[12] Inf
## ARIMA(0,0,3)(2,1,1)[12] 2125.264274
## ARIMA(0,0,0)(0,1,2)[12] 2437.395345
## ARIMA(1,0,0)(0,1,2)[12] 1954.192957
## ARIMA(2,0,0)(0,1,2)[12] 1937.290668
## ARIMA(3,0,0)(0,1,2)[12] 1933.592903
## ARIMA(4,0,0)(0,1,2)[12] 1932.720243
## ARIMA(0,0,1)(0,1,2)[12] 2249.260456
## ARIMA(1,0,1)(0,1,2)[12] 1929.409916
## ARIMA(2,0,1)(0,1,2)[12] 1930.353260
## ARIMA(3,0,1)(0,1,2)[12] 1933.572605
## ARIMA(0,0,2)(0,1,2)[12] 2157.535250
## ARIMA(1,0,2)(0,1,2)[12] 1929.232637
## ARIMA(2,0,2)(0,1,2)[12] 1932.181491
## ARIMA(0,0,3)(0,1,2)[12] 2115.301575
## ARIMA(1,0,3)(0,1,2)[12] 1931.317529
## ARIMA(0,0,4)(0,1,2)[12] 2089.541621
## ARIMA(0,0,0)(1,1,2)[12] 2445.765067
## ARIMA(1,0,0)(1,1,2)[12] 1952.629646
## ARIMA(2,0,0)(1,1,2)[12] 1935.094750
## ARIMA(3,0,0)(1,1,2)[12] 1933.956372
## ARIMA(0,0,1)(1,1,2)[12] 2253.554136
## ARIMA(1,0,1)(1,1,2)[12] 1927.883257
## ARIMA(2,0,1)(1,1,2)[12] 1931.413980
## ARIMA(0,0,2)(1,1,2)[12] 2160.445106
## ARIMA(1,0,2)(1,1,2)[12] Inf
## ARIMA(0,0,3)(1,1,2)[12] 2118.241676
## ARIMA(0,0,0)(2,1,2)[12] 2455.298789
## ARIMA(1,0,0)(2,1,2)[12] 1953.981451
## ARIMA(2,0,0)(2,1,2)[12] 1934.441987
## ARIMA(0,0,1)(2,1,2)[12] 2260.209505
## ARIMA(1,0,1)(2,1,2)[12] Inf
## ARIMA(0,0,2)(2,1,2)[12] 2169.280387
##
## --- Re-estimating best models without approximation ---
##
## ARIMA(2,0,1)(2,1,1)[12]+c 358.714599
tidy(fit)
## # A tibble: 18 × 8
## State Industry .model term estimate std.error statistic p.value
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 South Australia Hardware… arima… ar1 0.943 0.0209 45.1 7.01e-141
## 2 South Australia Hardware… arima… ma1 -0.286 0.0628 -4.56 7.35e- 6
## 3 South Australia Hardware… arima… sar1 -0.622 0.0565 -11.0 4.00e- 24
## 4 South Australia Hardware… arima… sar2 -0.266 0.0585 -4.54 7.90e- 6
## 5 South Australia Hardware… arima… cons… 0.0417 0.0169 2.47 1.40e- 2
## 6 South Australia Hardware… arima… ar1 1.09 0.166 6.55 2.30e- 10
## 7 South Australia Hardware… arima… ar2 -0.131 0.153 -0.858 3.91e- 1
## 8 South Australia Hardware… arima… ma1 -0.417 0.151 -2.75 6.30e- 3
## 9 South Australia Hardware… arima… sar1 -0.623 0.0564 -11.0 3.05e- 24
## 10 South Australia Hardware… arima… sar2 -0.262 0.0586 -4.47 1.08e- 5
## 11 South Australia Hardware… arima… cons… 0.0325 0.0137 2.37 1.83e- 2
## 12 South Australia Hardware… arima… ar1 1.13 0.149 7.60 3.24e- 13
## 13 South Australia Hardware… arima… ar2 -0.167 0.139 -1.20 2.29e- 1
## 14 South Australia Hardware… arima… ma1 -0.478 0.132 -3.62 3.39e- 4
## 15 South Australia Hardware… arima… sar1 0.00386 0.0764 0.0504 9.60e- 1
## 16 South Australia Hardware… arima… sar2 0.0234 0.0719 0.325 7.45e- 1
## 17 South Australia Hardware… arima… sma1 -0.797 0.0555 -14.4 1.99e- 36
## 18 South Australia Hardware… arima… cons… 0.0132 0.00262 5.02 8.53e- 7
glance(fit)
## # A tibble: 3 × 10
## State Industry .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
## 1 South Aust… Hardwar… arima… 0.191 -191. 393. 393. 416. <cpl> <cpl>
## 2 South Aust… Hardwar… arima… 0.192 -190. 394. 395. 421. <cpl> <cpl>
## 3 South Aust… Hardwar… arima… 0.167 -171. 358. 359. 388. <cpl> <cpl>
fit |>
forecast(h=24) |>
autoplot(train, level=NULL)
fit |>
forecast(h=24) |>
accuracy(myseries)
## # A tibble: 3 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 arima101… Sout… Hardwar… Test -6.24 9.08 8.01 -7.58 9.54 1.43 1.11 0.518
## 2 arima201… Sout… Hardwar… Test -6.13 8.91 7.87 -7.46 9.38 1.40 1.09 0.505
## 3 arima_au… Sout… Hardwar… Test -3.59 6.85 6.06 -4.41 7.12 1.08 0.840 0.445
ARIMA p,d,q-
myseries |>
gg_tsdisplay(difference(box_cox(Turnover, lambda)), plot_type = "partial")
myseries |>
model(ARIMA(box_cox(Turnover, lambda))) |>
tidy()
## # A tibble: 5 × 8
## State Industry .model term estimate std.error statistic p.value
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 South Australia Hardware,… ARIMA… ar1 0.972 0.0124 78.6 6.60e-257
## 2 South Australia Hardware,… ARIMA… ma1 -0.254 0.0496 -5.13 4.39e- 7
## 3 South Australia Hardware,… ARIMA… ma2 -0.143 0.0461 -3.09 2.14e- 3
## 4 South Australia Hardware,… ARIMA… sma1 -0.767 0.0326 -23.5 2.81e- 79
## 5 South Australia Hardware,… ARIMA… cons… 0.00889 0.00289 3.07 2.27e- 3
myseries |>
model(auto= ARIMA(box_cox(Turnover, lambda), trace=FALSE) ,
ets=ETS(box_cox(Turnover, lambda)~ trend("M"))
)|>
forecast(h=10, time=100) |>
autoplot(myseries)
When multiplicative trend in ETS, the plots are very similar.
myseries |>
model(auto= ARIMA(box_cox(Turnover, lambda), trace=FALSE) ,
ets=ETS(box_cox(Turnover, lambda)~ trend("M"))
)|>
accuracy()
## # A tibble: 2 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 South Aust… Hardwar… auto Trai… -0.0181 3.67 2.59 -0.250 5.49 0.421 0.427
## 2 South Aust… Hardwar… ets Trai… -0.284 4.13 2.95 -0.472 6.17 0.480 0.480
## # ℹ 1 more variable: ACF1 <dbl>
myseries |>
model(auto= ARIMA(box_cox(Turnover, lambda), trace=FALSE) ,
ets=ETS(box_cox(Turnover, lambda)~ trend("M"))
)|>
select(ets) |>
gg_tsresiduals()
Australian retail trade turnover. rdrr.io. https://rdrr.io/cran/tsibbledata/man/aus_retail.html
Chang, W. (2023, May 16). 13.13 Creating a QQ Plot | R Graphics Cookbook, 2nd edition. https://r-graphics.org/recipe-miscgraph-qq
Coder, R. (2021). Histogram with density in ggplot2. R CHARTS | a Collection of Charts and Graphs Made With the R Programming Language. https://r-charts.com/distribution/histogram-density-ggplot2/
Forecasting: Principles and Practice (3rd ed). (n.d.). https://otexts.com/fpp3/
Holtz, Y. (n.d.). The R Graph Gallery – Help and inspiration for R charts. The R Graph Gallery. https://r-graph-gallery.com/
Retail Trade, Australia, March 2023. (2023, May 9). Australian Bureau of Statistics. https://www.abs.gov.au/statistics/industry/retail-and-wholesale-trade/retail-trade-australia/latest-release
Robjhyndman. (n.d.). fpp3package/README.Rmd at master · robjhyndman/fpp3package. GitHub. https://github.com/robjhyndman/fpp3package/blob/master/README.Rmd
Roy, A. (2023). What is Sales Forecasting? CX Today. https://www.cxtoday.com/contact-centre/what-is-sales-forecasting/